#include "Calculus.h"

#define MIN(a, b) ((a) < (b) ? (a) : (b))
#define MAX(a, b) ((a) > (b) ? (a) : (b))
#define ACCURACY 1e-8
#define ITERATION 1000

using namespace Calculus;
using namespace LinearAlgebra;

// 显式构造函数
Polynomial::Polynomial(int o) : m_coefficients(std::vector<Real>(o, 0))
{
}

// 初始化列表构造
Polynomial::Polynomial(std::initializer_list<Real> l)
{
    for (auto it = l.begin(); it != l.end();)
        m_coefficients.push_back(*it++);
}

// 使用向量构造
Polynomial::Polynomial(std::vector<Real> &v) : m_coefficients(v)
{
}

// 拷贝构造
Polynomial::Polynomial(const Polynomial &poly) : m_coefficients(poly.m_coefficients)
{
}

// 多项式赋值
void Polynomial::operator=(const Polynomial &poly)
{
    m_coefficients = poly.m_coefficients;
}

// 加减法
Polynomial Polynomial::operator+(const Polynomial &poly) const
{
    int index = 0;
    int O1 = order();
    int O2 = poly.order();

    // 取最大和最小的阶数
    int Maxsize = MAX(O1, O2);
    int Minsize = MIN(O1, O2);

    // 新的多项式阶数是最大的那个
    Polynomial res(Maxsize);

    // 先计算两者都有的次数的项
    while (index < Minsize)
    {
        res.m_coefficients[index] = m_coefficients[index] + poly.m_coefficients[index];
        index++;
    }

    // 添加补充的次数的项
    while (index < O1)
    {
        res.m_coefficients[index] = m_coefficients[index];
        index++;
    }
    while (index < O2)
    {
        res.m_coefficients[index] = res.m_coefficients[index] + poly.m_coefficients[index];
        index++;
    }

    return res;
}
Polynomial Polynomial::operator-(const Polynomial &poly) const
{
    int index = 0;
    int O1 = order();
    int O2 = poly.order();

    // 取最大和最小的阶数
    int Maxsize = MAX(O1, O2);
    int Minsize = MIN(O1, O2);

    // 新的多项式阶数是最大的那个
    Polynomial res(Maxsize);

    // 先计算两者都有的次数的项
    while (index < Minsize)
    {
        res.m_coefficients[index] = m_coefficients[index] - poly.m_coefficients[index];
        index++;
    }

    // 添加补充的次数的项
    while (index < O1)
    {
        res.m_coefficients[index] = m_coefficients[index];
        index++;
    }
    while (index < O2)
    {
        res.m_coefficients[index] = res.m_coefficients[index] - poly.m_coefficients[index];
        index++;
    }

    return res;
}

// 乘常数
Polynomial Polynomial::operator*(Real k) const
{
    int O = order();
    Polynomial res(O);
    for (int i = 0; i < O; i++)
        res.m_coefficients[i] = m_coefficients[i] * k;

    return res;
}
Polynomial Polynomial::operator/(Real k) const
{
    assert(fabs(k) > ACCURACY);
    int O = order();
    Polynomial res(O);
    for (int i = 0; i < O; i++)
        res.m_coefficients[i] = m_coefficients[i] / k;

    return res;
}

// 多项式相乘
Polynomial Polynomial::operator*(const Polynomial &poly) const
{
    const int O1 = order();
    const int O2 = poly.order();
    Polynomial res(O1 + O2 - 1);
    for (int i = 0; i < O1; i++)
        for (int j = 0; j < O2; j++)
            res.m_coefficients[i + j] += m_coefficients[i] * poly.m_coefficients[j];

    return res;
}

// 带余除法，获得商多项式
Polynomial Polynomial::operator/(const Polynomial &poly) const
{
    int O1 = order();
    int O2 = poly.order();

    // 阶数不够，返回 0 多项式
    if (O1 < O2)
        return Polynomial(0);

    Polynomial res(O1 - O2 + 1);

    // 最高阶系数
    res.m_coefficients[O1 - O2] = m_coefficients[O1 - 1] / poly.m_coefficients[O2 - 1];

    // 获得余式
    Polynomial tmp = *this - res * poly;

    // 减去最高项后需要降阶，然后递归计算商
    tmp.reduce();
    tmp = tmp / poly;

    // 商式赋值
    for (int i = 0; i < tmp.order(); ++i)
        res.m_coefficients[i] = tmp.m_coefficients[i];

    return res;
}

// 多项式求值
Real Polynomial::operator()(Real x) const
{
    int O = order();
    Real y = 0, xi = 1;
    for (int i = 0; i < O; i++)
    {
        y += m_coefficients[i] * xi;
        xi *= x;
    }
    return y;
}

// 推入高阶项
void Polynomial::push(Real k)
{
    m_coefficients.push_back(k);
}

// 多项式降阶，将系数为 0 的高阶项去除
void Polynomial::reduce()
{
    // 只要最后一项非常小，就弹出最后一个元素，除非没有元素
    while (order() > 0 && fabs(m_coefficients.back()) < ACCURACY)
        m_coefficients.pop_back();
}

// 定积分
Real Polynomial::intergral(Real a, Real b) const
{
    int O = order();
    Real I = 0;
    for (int i = 0; i < O; i++)
        I += m_coefficients[i] * (pow(b, i + 1) - pow(a, i + 1)) / (i + 1);

    return I;
}

// 不定积分
Polynomial Polynomial::intergral() const
{
    Polynomial poly;
    poly.push(0);
    for (int i = 0; i < order(); i++)
        poly.push(coeff(i) / (i + 1));

    return poly;
}

// 导函数
Polynomial Polynomial::derivative() const
{
    Polynomial poly;
    for (int i = 1; i < order(); i++)
        poly.push(coeff(i) * i);

    return poly;
}

// 获取多项式的阶
int Polynomial::order() const
{
    return m_coefficients.size();
}

// 获得对应次数的系数
Real Polynomial::coeff(int n) const
{
    assert(n < m_coefficients.size());
    return m_coefficients[n];
}

// 多项式的拉普拉斯阵
Matrix Polynomial::Laplace() const
{
    int n = order() - 1;
    Matrix A(n, n);
    for (int i = 0; i < n; i++)
    {
        A(0, n - i - 1) = -m_coefficients[i];
        if (i < n - 1)
            A(i + 1, i) = 1;
    }
    return A;
}

// 没有区间端点，只有一个多项式
PiecewisePolynomial::PiecewisePolynomial(Polynomial poly) : m_poly({poly})
{
}

// 拷贝构造
PiecewisePolynomial::PiecewisePolynomial(const PiecewisePolynomial &ppoly) : m_interval(ppoly.m_interval), m_poly(ppoly.m_poly)
{
}

// 只设置区间，每段区间上是零
PiecewisePolynomial::PiecewisePolynomial(std::deque<Real> &interval) : m_interval(interval), m_poly(std::deque<Polynomial>(interval.size() + 1))
{
}

// 设置区间和多项式
PiecewisePolynomial::PiecewisePolynomial(std::deque<Real> &interval, std::deque<Polynomial> &poly) : m_interval(interval), m_poly(poly)
{
    // 检查赋值是否有效
    assert(interval.size() + 1 == poly.size());
}

// 多项式赋值
void PiecewisePolynomial::operator=(const PiecewisePolynomial &poly)
{
    m_interval = poly.m_interval;
    m_poly = poly.m_poly;
}

// 多项式加减法
PiecewisePolynomial PiecewisePolynomial::operator+(const PiecewisePolynomial &poly) const
{
    // 必须具有相同的区间端点
    assert(m_interval.size() == poly.m_interval.size());
    for (int i = 0; i < m_interval.size(); ++i)
        assert(fabs(m_interval[i] - poly.m_interval[i]) < ACCURACY);

    PiecewisePolynomial res = *this;
    for (int i = 0; i < m_poly.size(); ++i)
        res.m_poly[i] = res.m_poly[i] + poly.m_poly[i];

    return res;
}
PiecewisePolynomial PiecewisePolynomial::operator-(const PiecewisePolynomial &poly) const
{
    // 必须具有相同的区间端点
    assert(m_interval.size() == poly.m_interval.size());
    for (int i = 0; i < m_interval.size(); ++i)
        assert(fabs(m_interval[i] - poly.m_interval[i]) < ACCURACY);

    PiecewisePolynomial res = *this;
    for (int i = 0; i < m_poly.size(); ++i)
        res.m_poly[i] = res.m_poly[i] - poly.m_poly[i];

    return res;
}
PiecewisePolynomial PiecewisePolynomial::operator+(const Polynomial &poly) const
{
    PiecewisePolynomial res = *this;
    for (int i = 0; i < m_poly.size(); ++i)
        res.m_poly[i] = res.m_poly[i] + poly;

    return res;
}
PiecewisePolynomial PiecewisePolynomial::operator-(const Polynomial &poly) const
{
    PiecewisePolynomial res = *this;
    for (int i = 0; i < m_poly.size(); ++i)
        res.m_poly[i] = res.m_poly[i] - poly;

    return res;
}

// 多项式乘常数
PiecewisePolynomial PiecewisePolynomial::operator*(Real k) const
{
    PiecewisePolynomial res = *this;
    for (int i = 0; i < m_poly.size(); ++i)
        res.m_poly[i] = res.m_poly[i] * k;

    return res;
}
PiecewisePolynomial PiecewisePolynomial::operator/(Real k) const
{
    PiecewisePolynomial res = *this;
    for (int i = 0; i < m_poly.size(); ++i)
        res.m_poly[i] = res.m_poly[i] / k;

    return res;
}

// 多项式相乘，分段相乘
PiecewisePolynomial PiecewisePolynomial::operator*(const PiecewisePolynomial &poly) const
{
    // 否则必须具有相同的区间端点
    assert(m_interval.size() == poly.m_interval.size());
    for (int i = 0; i < m_interval.size(); ++i)
        assert(fabs(m_interval[i] - poly.m_interval[i]) < ACCURACY);

    PiecewisePolynomial res = *this;
    for (int i = 0; i < m_poly.size(); ++i)
        res.m_poly[i] = res.m_poly[i] * poly.m_poly[i];

    return res;
}
PiecewisePolynomial PiecewisePolynomial::operator*(const Polynomial &poly) const
{
    PiecewisePolynomial res = *this;
    for (int i = 0; i < m_poly.size(); ++i)
        res.m_poly[i] = res.m_poly[i] * poly;

    return res;
}

// 多项式除法
PiecewisePolynomial PiecewisePolynomial::operator/(const Polynomial &poly) const
{
    PiecewisePolynomial res = *this;
    for (int i = 0; i < m_poly.size(); ++i)
        res.m_poly[i] = res.m_poly[i] / poly;

    return res;
}

// 多项式求值
Real PiecewisePolynomial::operator()(Real x) const
{
    return poly(x)(x);
}

// 在右端添加端点，增加一个零多项式
void PiecewisePolynomial::push_back(Real t)
{
    // 如果区间端点不是空，就检查是否大于最大的端点，如果重合或小于，就什么也不做
    if (!m_interval.empty() && t < m_interval.back() + ACCURACY)
        return;

    m_interval.push_back(t);
    m_poly.push_back(Polynomial());
}

// 在左端添加端点，增加一个零多项式
void PiecewisePolynomial::push_front(Real t)
{
    if (!m_interval.empty() && t > m_interval.front() - ACCURACY)
        return;

    m_interval.push_front(t);
    m_poly.push_front(Polynomial());
}

// 获得指定索引的多项式
Polynomial PiecewisePolynomial::operator[](int i) const
{
    assert(i < m_poly.size());
    return m_poly[i];
}

// 获得指定索引的多项式引用
Polynomial &PiecewisePolynomial::operator[](int i)
{
    assert(i < m_poly.size());
    return m_poly[i];
}

// 获得传入点对应的多项式
Polynomial PiecewisePolynomial::poly(Real t) const
{
    int index = 0;

    // 找到 x 所属的区间范围的多项式
    for (int i = 0; i < m_interval.size(); ++i)
    {
        if (t > m_interval[i])
            index++;
    }

    return m_poly[index];
}

// 获得端点数
int PiecewisePolynomial::size() const
{
    return m_interval.size();
}

// 获得指定索引的阶数
int PiecewisePolynomial::order(int i) const
{
    assert(i < m_poly.size());
    return m_poly[i].order();
}

// 不定积分
PiecewisePolynomial PiecewisePolynomial::intergral() const
{
    PiecewisePolynomial ppoly;
    ppoly.m_interval = m_interval;

    // 注意默认构造的分段多项式有一个零多项式，需要清除
    ppoly.m_poly.clear();
    for (int i = 0; i < m_poly.size(); i++)
        ppoly.m_poly.push_back(m_poly[i].intergral());

    return ppoly;
}

// 导函数
PiecewisePolynomial PiecewisePolynomial::derivative() const
{
    PiecewisePolynomial ppoly;
    ppoly.m_interval = m_interval;

    // 注意默认构造的分段多项式有一个零多项式，需要清除
    ppoly.m_poly.clear();
    for (int i = 0; i < m_poly.size(); i++)
        ppoly.m_poly.push_back(m_poly[i].derivative());

    return ppoly;
}

// 检查端点处的连续性
int PiecewisePolynomial::continuity(int i) const
{
    assert(i < m_interval.size());

    // 左右两个多项式
    Polynomial left = m_poly[i];
    Polynomial right = m_poly[i + 1];

    int min = MIN(left.order(), right.order());
    Real t = m_interval[i];

    // 标记连续性，-1 表示不连续
    int c = -1;
    for (int j = 0; j < min; j++)
    {
        if (fabs(left(t) - right(t)) < ACCURACY)
            c++;
        left = left.derivative();
        right = right.derivative();
    }

    return c;
}

// 整体连续性
int PiecewisePolynomial::continuity() const
{
    // 不允许检查纯多项式的连续性，其连续性直接是无穷大
    assert(m_interval.size() > 0);

    int c = continuity(0);
    for (int i = 1; i < m_interval.size(); i++)
        c = MIN(c, continuity(i));
    return c;
}

// 合并两个分段多项式
PiecewisePolynomial Calculus::merge(const PiecewisePolynomial &p1, const PiecewisePolynomial &p2)
{
    // 两个都是纯多项式，没有端点，直接相加
    if (p1.size() == 0 && p2.size() == 0)
        return p1 + p2;

    // 首先合并区间端点
    PiecewisePolynomial res;
    int n1 = p1.size();
    int n2 = p2.size();

    // 储存所有端点，然后排序
    Vector x(n1 + n2);
    for (int i = 0; i < n1; ++i)
        x[i] = p1.m_interval[i];
    for (int i = 0; i < n2; ++i)
        x[i + n1] = p2.m_interval[i];
    x.sort();

    // 推入这些端点
    for (int i = 0; i < n1 + n2; ++i)
        res.push_back(x[i]);

    // 开始计算加法，我们保存每个区间的中点，左右添加两个点
    std::deque<Real> mid;
    for (int i = 0; i < res.size() - 1; i++)
        mid.push_back((res.m_interval[i + 1] + res.m_interval[i]) / 2);
    mid.push_front(res.m_interval.front() - 1);
    mid.push_back(res.m_interval.back() + 1);

    // 获取对应的多项式相加
    for (int i = 0; i < mid.size(); i++)
        res.m_poly[i] = p1.poly(mid[i]) + p2.poly(mid[i]);

    return res;
}

// 多元向量微分，返回 Jacobi 矩阵
Matrix Jacobian::operator()(FuncXX f, const Vector &x)
{
    Vector y = f(x);
    int m = y.size();
    int n = x.size();
    Real h = 1e-8;

    Matrix d(m, n);
    for (int i = 0; i < n; i++)
    {
        Vector y1 = x, y2 = x;
        y1[i] += h, y2[i] -= h;

        // 中心差分向量
        Vector dp = (f(y1) - f(y2)) / (2 * h);
        for (int j = 0; j < m; j++)
            d(j, i) = dp[j];
    }
    return d;
}

// 多元函数的 Jacobi 矩阵：由对每个分量的偏导构成的矩阵
Matrix Jacobian::operator()(FuncX1 f, const Vector &x)
{
    int n = x.size();
    Real h = 1e-8;

    // 1行矩阵
    Matrix d(1, n);
    for (int i = 0; i < n; i++)
    {
        Vector y1 = x, y2 = x;
        y1[i] += h, y2[i] -= h;

        // 中心差分计算导数
        d(0, i) = (f(y1) - f(y2)) / (2 * h);
    }
    return d;
}

// 一元向量函数的导向量：每个分量的导数构成的向量
Vector Jacobian::operator()(Func1X f, Real x)
{
    Real h = 1e-8;

    // 直接计算向量作差
    Vector y1 = f(x + h);
    Vector y2 = f(x - h);

    return (y1 - y2) / (2 * h);
}

// 一元函数的导数
Real Jacobian::operator()(Func11 f, Real x)
{
    Real h = 1e-8;
    return (f(x + h) - f(x - h)) / (2 * h);
}

// 初始化权
Integrator::Integrator(int N) : m_weight(Vector(N + 1))
{
}

// 积分器分割数
int Integrator::order() const
{
    return m_weight.size() - 1;
}

// 初始化权重
NewtonCotes::NewtonCotes(int N) : Integrator(N)
{
    // 使用插值多项式积分计算权
    Vector x(N + 1);
    for (int i = 0; i < N + 1; i++)
        x[i] = i;

    for (int i = 0; i < N + 1; i++)
    {
        Vector f(N + 1);
        f[i] = 1;
        m_weight[i] = interpolate(x, f).intergral(0, N);
    }
}

// 牛顿-柯茨公式
Real NewtonCotes::operator()(Func11 f, Real a, Real b)
{
    int N = m_weight.size() - 1;

    // 分割宽度
    Real h = (b - a) / N;

    // 节点值
    Vector y(N + 1);
    for (int i = 0; i < N + 1; i++)
        y[i] = f(h * i + a);

    // 节点加权
    return dot(m_weight, y) * h;
}

// 初始化权节点和权重
GaussIntegrator::GaussIntegrator(int N) : Integrator(N), m_node(Vector(N + 1))
{
    // 获得 Gauss-Legendre 多项式
    Polynomial GL = Legendre(N + 1);

    // 使用牛顿法迭代求根
    FSolve FS;

    // 依次计算 N+1 个根，从 0 开始迭代，通常会得到从小到大的根
    for (int i = 0; i < N + 1; i++)
    {
        // 构造函数
        Func11 F = [=](Real x)
        {
            return GL(x);
        };

        // 迭代求根，从 1 开始，因为 Legendre 多项式的根不在 -1,1 之外，防止初始值与根重合，或者出现导数为零的情况
        m_node[i] = FS(F, 1);

        // 除以因式 x - c[i] 进行降阶
        GL = GL / Polynomial({-m_node[i], 1});
    }

    // 将根按照升序排列
    m_node.sort();

    // 初始化权，在 -1,1 的正交多项式区间上对初等拉格朗日函数积分
    for (int i = 0; i < N + 1; i++)
        m_weight[i] = Lagrange(m_node, i).intergral(-1, 1);
}

// 高斯积分公式
Real GaussIntegrator::operator()(Func11 f, Real a, Real b)
{
    int N = m_weight.size() - 1;

    // 节点值，需要区间从 -1,1 变换到 a,b
    Vector y(N + 1);
    for (int i = 0; i < N + 1; i++)
        y[i] = (b - a) / 2 * f((b - a) / 2 * m_node[i] + (b + a) / 2);

    // 节点加权
    return dot(m_weight, y);
}

// 一元函数求根
Real FSolve::operator()(Func11 f, Real x)
{
    // Jacobian 算子
    Jacobian j;

    Real dx = 0;
    for (int i = 0; i < ITERATION; i++)
    {
        Real fx = f(x);

        // 获得导数
        Real J = j(f, x);

        // 迭代推进 x = x - f / f'
        dx = fx / J;
        x = x - dx;

        // 当变化量足够小或者是根就退出
        if (fabs(dx) < 1e-8 || fabs(fx) < 1e-8)
            break;
    }

    return x;
}

// 多元向量函数牛顿法求根
void FSolve::operator()(FuncXX f, Vector &x)
{
    // Jacobian 算子
    Jacobian j;

    // 线性求解器
    LinearSolver solver;

    // 创建容器
    Vector dx(x.size());
    for (int i = 0; i < ITERATION; i++)
    {
        // 获得 Jacobi 矩阵
        Matrix J = j(f, x);
        if (J.norm() < ACCURACY)
            break;

        // JT*Jdx = -JT*f
        Vector fx = J.transpose() * f(x);
        Matrix JTJ = J.transpose() * J;

        // 求解最小二乘问题
        dx = solver.solve(JTJ, fx, true);
        x = x - dx;

        // 当变化量足够小或者是根就退出
        if (dx.norm() < ACCURACY || fx.norm() < ACCURACY)
            break;
    }
}

// 阶乘
Real Calculus::factorial(int n)
{
    Real i = 1;
    for (int k = 0; k < n; k++)
        i *= k + 1;

    return i;
}

// 组合数
Real Calculus::combine(int n, int i)
{
    return factorial(n) / (factorial(i) * factorial(n - i));
}

// 多项式整数幂
Polynomial Calculus::power(const Polynomial &poly, int n)
{
    Polynomial res({1});
    for (int i = 0; i < n; i++)
        res = res * poly;
    return res;
}

// 产生节点多项式
Polynomial Calculus::nodes(Vector &c)
{
    Polynomial poly({1});
    for (int i = 0; i < c.size(); i++)
        poly = poly * Polynomial({-c[i], 1});

    return poly;
}

// 已知节点 c0 c1 ... cn ，获得 cj 对应的初等多项式
Polynomial Calculus::Lagrange(Vector &c, int j)
{
    Polynomial poly({1});
    for (int k = 0; k < c.size(); k++)
        if (k != j)
            poly = poly * Polynomial({-c[k], 1}) / (c[j] - c[k]);

    return poly;
}

// 多项式插值
Polynomial Calculus::interpolate(Vector &x, Vector &f)
{
    int n = f.size();
    std::vector<std::vector<Real>> table; // 差分表

    // 计算差分表
    for (int i = 0; i < n; i++)
    {
        table.push_back({f[i]});
        for (int j = 1; j <= i; j++)
            table[i].push_back((table[i][j - 1] - table[i - 1][j - 1]) / (x[i] - x[i - j]));
    }

    // 由差分表求取插值多项式
    Polynomial poly(n);
    for (int i = 0; i < table.size(); i++)
    {
        Polynomial tmp({table[i][i]});

        // 多项式连乘
        for (int j = i - 1; j >= 0; j--)
            tmp = tmp * Polynomial({-x[j], 1});

        poly = poly + tmp;
    }

    return poly;
}

// Gauss-Legendre 多项式
Polynomial Calculus::Legendre(int s)
{
    Polynomial F0({1});
    if (s > 0)
    {
        Polynomial x({0, 1});
        Polynomial F1 = x;
        Polynomial F2 = x;
        for (int i = 1; i < s; i++)
        {
            F2 = x * F1 * (2 * i + 1) / (i + 1) - F0 * i / (i + 1);
            F0 = F1;
            F1 = F2;
        }
        F0 = F2;
    }
    return F0;
}

// 获得 Bernstein 基函数 Cni * t^i * (1-t)^(n-i)
Polynomial Calculus::Bernstein(int n, int i)
{
    assert(n >= i);
    return power(Polynomial({0, 1}), i) * power(Polynomial({1, -1}), n - i) * combine(n, i);
}

// 获得 B 样条基函数
PiecewisePolynomial Calculus::BasicBspline(std::deque<Real> &knot)
{
    // 至少要两个节点才能计算
    assert(knot.size() > 1);

    // 递归到最底层，恰好有两个节点
    if (knot.size() == 2)
    {
        // 如果两个节点重合，返回纯零多项式
        if (knot[0] == knot[1])
            return PiecewisePolynomial();

        // 如果不重合，设置节点之间取常数 1
        PiecewisePolynomial Bik(knot);
        Bik[1] = Polynomial({1});
        return Bik;
    }

    // 减少一个节点得到第一阶的两个样条
    std::deque<Real> knot1 = knot;
    std::deque<Real> knot2 = knot;
    knot1.pop_back();
    knot2.pop_front();

    // 计算低阶样条
    PiecewisePolynomial Bik_1 = BasicBspline(knot1);
    PiecewisePolynomial Bi_1k_1 = BasicBspline(knot2);

    // 递推公式 B_i,k = (t-t_i)/(t_i+k-1-t_i)B_i,k-1 + (t_i+k - t)/(t_i+k - t_i+1)B_i+1,k-1
    Real t_i = knot[0], t_i1 = knot[1];
    Real t_ik_1 = knot[knot.size() - 2], t_ik = knot[knot.size() - 1];

    // 如果分母为零，则规定系数为零
    Real coeff1 = (fabs(t_ik_1 - t_i) < ACCURACY) ? 0 : (1 / (t_ik_1 - t_i));
    Real coeff2 = (fabs(t_ik - t_i1) < ACCURACY) ? 0 : (1 / (t_ik - t_i1));

    // 最后合并两个多项式
    return merge(Bik_1 * Polynomial({-t_i, 1}) * coeff1, Bi_1k_1 * Polynomial({t_ik, -1}) * coeff2);
}